#!/usr/bin/perl
# 2013/09/10 Extract barcode and add to ID
# 2012/11/15 Split sequence reads by detecint barcode using BLAST

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use File::Basename;
use Carp qw(croak);
$| = 0;

if(@ARGV != 2 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [original fastq]\n";
}

#---------------------------------------
# parameter
#---------------------------------------
my %opt;
getopts("i:", \%opt);
my $FILE_original = $opt{i};

my $LINKER_LENGTH = 34;	# Both ChIA-PET and HiD are same
my ($name, $dir, $ext) = &fileparse($FILE_original, '\..*');
my $FILE_table = $dir . $name . '_table.txt';
my $FILE_fasta = $dir . $name . '.fasta';
my $FILE_blast = $dir . $name . '_blast_result.txt';
my $FILE_fastq1 = $dir . $name . '_1.fastq';
my $FILE_fastq2 = $dir . $name . '_2.fastq';


#---------------------------------------
# flow chart of program
#---------------------------------------
# 1. making fasta file
# 2. BLAST to search linker seqeunce
# 3. split file based on the linker position


#---------------------------------------
# making fasta file
#---------------------------------------
{
	&Progress("Output fasta file");
	my $fh_original = IO::File->new($FILE_original) or die "cannot open $FILE_original: $!";
	my $fh_table = IO::File->new($FILE_table, 'w') or die "cannot write $FILE_table: $!";
	my $fh_fasta = IO::File->new($FILE_fasta, 'w') or die "cannot write $FILE_fasta: $!";

	my $line = 0;
	my $id;
	my $seq;
	my $quality;
	my $sequenceNum = 0;
	while($_ = $fh_original->getline()){
		s/\r?\n//;
		$line++;
		if($line % 4 == 1){
			$sequenceNum++;
			$id = $sequenceNum;
			if($sequenceNum % 5000000 == 0){
				&Progress("$sequenceNum finished...");
			}
		}elsif($line % 4 == 2){
			$seq = $_;
		}elsif($line % 4 == 0){
			$quality = $_;

			# output
			$fh_fasta->print(">$id\n$seq\n");
			$fh_table->print("$id\t$seq\t$quality\n");
		}
	}
	$fh_original->close();
	$fh_table->close();
	$fh_fasta->close();
}



#---------------------------------------
# BLAST
#---------------------------------------
&Progress("Blast start");
my $blast_CMD =<<"EOF";
blastn -query $FILE_fasta -strand 'plus' -db HiD_linker -task blastn -out $FILE_blast -outfmt "6 qseqid sseqid qstart qend sstart send" -max_target_seqs 1;
EOF

### example output
#ID_0	ChIA-PET_linkerBB	24	57	1	34
#ID_1	ChIA-PET_linkerBB	24	57	1	34
#ID_2	ChIA-PET_linkerBB	24	57	1	34
#ID_3	ChIA-PET_linkerBB	26	57	3	34
#ID_4	ChIA-PET_linkerBB	24	57	1	34
#ID_5	ChIA-PET_linkerBB	24	57	1	34
#ID_6	ChIA-PET_linkerBB	24	57	1	34
#ID_7	ChIA-PET_linkerBB	24	57	1	34
#ID_8	ChIA-PET_linkerBB	24	57	1	34
#ID_9	ChIA-PET_linkerBB	24	55	1	32

system($blast_CMD) == 0 or die "BLAST failed; $!";



#---------------------------------------
# split file
#---------------------------------------
{
	&Progress("Split file start");
	my $fh_blast = IO::File->new($FILE_blast) or die "cannot open $FILE_blast: $!";
	my $fh_table = IO::File->new($FILE_table) or die "cannot open $FILE_table: $!";
	my $fh_seq1 = IO::File->new($FILE_fastq1, 'w') or die "cannot write $FILE_fastq1: $!";
	my $fh_seq2 = IO::File->new($FILE_fastq2, 'w') or die "cannot write $FILE_fastq2: $!";

	while(my $blast = $fh_blast->getline()){
		my $table = $fh_table->getline();
		$blast =~ s/\r?\n//;
		$table =~ s/\r?\n//;

		my ($bid, $linker, $qstart, $qend, $sstart, $send) = split /\t/, $blast;
		my ($tid, $seq, $quality) = split /\t/, $table;

		### sometimes blast cannot find linker sequnece and skipped those read
		### then skip treatment of those read and move to next
		while($tid < $bid){
			$table = $fh_table->getline();
			($tid, $seq, $quality) = split /\t/, $table;
		}

		if($bid != $tid){
			die ("blast and table id not matched ! at \n$blast\n$table\n");
		}

		if($bid % 5000000 == 0){
			&Progress("$bid finished...");
		}

		# check barcode sequence
		my $barcode = 'NA';
		my $barcode_target_seq = substr($seq, 84);	# at least read should have 25bp + 34bp + 25bp = 84bp in the left
		if($barcode_target_seq =~ m/GAC([ATGC]{5})GCA/){
			$barcode = $1;
		}

		# combin id
		my $ID = '@' . $bid . ':' . $barcode;

		# adjusting the alingment position in case of partial alignment
		if($sstart != 1){
			$qstart = $qstart - ($sstart - 1);
		}
		if($send != $LINKER_LENGTH){
			$qend = $qend + $LINKER_LENGTH - $send;
		}

		# start begins from 0
		$qstart -= 1;
		$qend -= 1;

		# remove sequence if remaing length truncated sequence become less than 20bp
		if(($qstart < 20) or ($qend  > 80)){
			next;
		}

		### frist half part
		my ($seq_1, $seq_2);
		my ($quality_1, $quality_2);
		{
			my $start = $qstart - 25;
			if($start < 0){
				$start = 0;
			}
			$seq_1 = substr($seq, $start, 25);
			$quality_1 = substr($quality, $start, 25);
			$fh_seq1->printf("%s\n%s\n+\n%s\n", $ID, $seq_1, $quality_1);
		}
		{
			my $start = $qend + 1;
			$seq_2 = substr($seq, $start, 25);
			$quality_2 = substr($quality, $start, 25);
			$fh_seq2->printf("%s\n%s\n+\n%s\n", $ID, $seq_2, $quality_2);
		}

	}
	$fh_blast->close();
	$fh_table->close();
	$fh_seq1->close();
	$fh_seq2->close();
}


sub Progress{
	my ($message) = @_;
	my $time = &time();
	warn "[$time] $message\n";
}


sub time{
	my ($sec, $min, $hour, $mday, $mon, $year) = localtime(time);
	my $time = sprintf ("%d/%0d %02d:%02d",$mon+1, $mday, $hour, $min);
	return $time;
}